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Abstract. The energy spectrum of cosmic rays in the range E ~ 10^^ eV 
to 6 X 10 -"^^ eV has been studied using the air Cherenkov hght detectors of 
the Yakutsk array. The total flux of photons produced by relativistic electrons 
(including positrons as well, hereafter) of extensive air showers in the atmosphere 
is used as the energy estimator of the primary particle initiating a shower. The 
resultant differential flux of cosmic rays exhibits, in accordance with previous 
measurements, a knee and ankle features at energies 3 x lO-*^^ and ~ lO-*^^ eV, 
respectively. A comparison of observational data with simulations is made in the 
knee and ankle regions in order to choose the models of galactic and extragalactic 
components of cosmic rays which describe better the energy spectrum measured. 
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1. Introduction 

The Cherenkov hght emitted in the atmosphere by relativistic electrons of extensive 
air showers (EASs) of cosmic rays (CRs) carries important information about the 
shower development and the primary CR particles. Well-known application of the 
air Cherenkov light technique is the 7-ray astronomy. But in this article we aim at 
another application of the technique - namely to measure Cherenkov photons from 
EAS initiated by CR of energy above 10^^ eV. These giant showers supply plenty of 
light so that one can detect it with unarmed photomultiplier tubes (PMT) triggered 
within pulse duration. 

Since the first observation by Galbraith and Jelley [1 and a systematic 
measurement of air Cherenkov light properties in the Pamir experiment [2 , a 
number of EAS arrays have been equipped with Cherenkov light detectors. It seems 
that the most durable and plentiful in Cherenkov light data is the Yakutsk array 
experiment [3, 4 . The total flux of light is used to estimate the primary energy in a 
model independent manner and the radial distribution of the light intensity at ground 
level is used to infer the position of shower maximum, Xmax^ in the atmosphere [3 . 

In this article we focus on the experimental data obtained in Yakutsk with 
Cherenkov light detectors aiming at the CR energy spectrum. The paper is structured 
as follows. In Section |2] we outline the general characteristics of the Yakutsk array 
experiment, while in Sections [3] and [4] particular properties of the Cherenkov light 
detectors are given: detector design and calibration of the signal. In Section [5] the 
measurement and monitoring of the atmospheric extinction of light is described. The 
measured lateral distribution of the Cherenkov light intensity at the observation level 
is given in Section [6] The method used to estimate the energy of the particle initiating 
the EAS is described in Section [7[ Resultant energy spectrum of CRs is discussed in 
Section |8] Our conclusions are set out in Section [9j In two Appendices additional 
material is given essential for the subjects considered. 

2. The Yakutsk array 

The Yakutsk array is located in Oktyomtsy near Yakutsk, Russia (61. 7^ A^, 129. 4^£^), 
100 m above sea level (1020 g/cm^). At present it consists of 58 ground-based and 
6 underground scintillation detector stations to measure charged particles (electrons 
and muons) and 48 detectors - PMTs in shuttered housing to observe the atmospheric 
Cherenkov light. During more than 30 years of lifetime the Yakutsk array has been 
re-configured several times, the total area covered by detectors was maximal about 
1990 {Seff ~ 17 km^), now it is Seff ~ 10 km^. In the central part of the array there 
is a denser domain with 100-250 m detector spacing. During the whole observation 
period approximately 10^ showers of the primary energy above 10^^ eV are detected; 
the three highest energy events selected with axes within the array area and zenith 
angle < 60^ have an energy E > 10^^ eV. 

The actual detector arrangement of the array is shown in Figure [l] Charged 
particle detectors of 2 m^ area are built in stations in couples; the Cherenkov light 
detectors - PMTs of 176 cm^ and 3 x 176 cm^ acceptance area, forms the medium, Ci 
(~ 500 m spacing), and the autonomous, C2 (50 to 200 m spacing), subsets. The latter 
was added in 1995 with the aim to study air showers in the energy range 10^^ — 10^^ 
eV via the Cherenkov light measurements [5 . 
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Figure 1. The detector arrangement of the Yakutsk array. Charged particle 
detectors (open circles), Cherenkov light detectors of the Ci subset (filled circles) 
and the C2 subset (filled triangles), and the muon detectors (squares) are shown. 
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Figure 2. Local area network of the Yakutsk array. 



All detectors/controllers and data processing units of the array are connected into 
the data handling network shown in figure [2] 



3. Air Cherenkov light detector design 

Charged particle detectors of the array were described in [3^. Here we will detail a 
Cherenkov light detector unit. It consists of a vertically mounted PMT (FEU-49B, 
15 cm diameter) with amplifier in a metal container blackened inside [6 . An upper 
hole provides < 55^ aperture (figure [3|. To protect the photocathode from sunlight 
the motorized light-proof lid is set. At night, all lids of the array can be commanded 
remotely to open. PMTs and amplifiers are powered around-the-clock to guarantee 
a stability of performance. When the lid is open, a fan blows with warm air to keep 
snow and dust out of the photocathode surface. 

There is a variant of detector with three PMTs in a housing which can operate 
independently or in summation of signal in order to increase acceptance area at the 
shower periphery. In addition, dedicated detectors were used to measure the shape 
and width of the Cherenkov signal from the shower. As an example, the pulse shape of 
the Cherenkov signal at the shower periphery is shown in figure [4j while the halfwidth 
of the signal as a function of radial distance is given in figure [5] 

The spectral sensitivity of the PMT used is shown by the dashed line in 
figure [g] together with air Cherenkov light spectrum (solid line) and the atmospheric 



Measuring EAS with Cherenkov light detectors 



4 




Figure 3. The air Cherenkov hght detector. 



Figure 4. Time (horizontal axis, ns) dependence of the Cherenkov hght signal 
(arbitrary units, vertical axis). 

transmission curve. 

Two kind of triggers were used to select the showers from the background at the 
Yakutsk array: produced by the scintihators and by Cherenkov hght detectors. In 
the first case a coincidence signal (if the particle density, p, is greater than 0.5 m~^ 
in two scintillators of each station within 2 /is) passes on to the central controller. 
Trigger-500 is then produced in the case of a coincident signal (in 40 /is) from three or 
more stations with ~ 500 m spacing (> sixfold coincidence). Similarly, trigger-1000 is 
produced by ~ 1 km spacing stations. After 1992 when 18 new stations were added, 
the array area is increased from 2.5 km^ to 7.2 km^ where trigger-500 operates. That 
is why we can deal with EAS in the energy range from 3 x 10^^ eV to 3 x 10^^ eV 
using the same trigger. 

The Ci subset of Cherenkov light detectors has no hardware trigger of its own. 
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Figure 5. The half width of the Cherenkov hght pulse as a function of the shower 
core distance. Experiment: points 3 ; theory: solid line with la errors (dashed 
lines) E]. 
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Figure 6. The spectral sensitivity of the PMT (dashed line). Air transmission of 
a cloudless atmosphere (dotted line) and Cherenkov emission spectrum (full line). 



Instead, the scintillator trigger controls the detection of Cherenkov signal in this 
subarray. The amplitude of PMT signal is recorded within 5 /is after the coincident 
signal in two scintillators of the host station. In order to exclude detection of 
an accidental Cherenkov signal when there is no trigger-500/1000, a double trigger 
condition is in use in data analysis: 

• three or more scintillator stations have detected the charged particles {p > 0.5 



m 



0; 



• the Cherenkov light intensities in three or more PMTs are greater than 2.4 x 
10^/4.5 X 10^ m~^ depending on the acceptance of PMT. 
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The thresholds have been chosen to keep the signal/noise ratio greater than 3. 

The second case is C2 subset which has an independent Cherenkov hght trigger 
formed by three or more PMTs having detected hght intensities above a given 
threshold within 10 /is. The signal integration time of the individual PMT is 0.5 

/iS. 



4. Detector calibration 



In order to find the absolute response of the array to Cherenkov light we have 
considered the sensitivity of all the PMTs in the detectors and have measured the 
conversion factors from light intensity into ADC value. The Cherenkov light emitted 
by relativistic muons in a transparent medium with known refraction index (optical 
radiator) was used in the detector calibration, as was suggested in [8 . In this paper we 
are keeping along the description given in [6l |9] . Distilled water and plexiglas are used 
to measure the output signal in a geometry shown in figure [?[ A telescope composed 
of two PMTs above and below the radiator, spaced 1 m from the detector, bottom 
one under the lead shield (10 cm), selects high-energy muons from background CRs. 
On the photocathode of the detector PMT is placed (above a thin layer of glycerin) 
a plexiglas disk of 5.5 cm thickness and 15 cm diameter with polished side surface. 
As another version of radiator, distilled water was used without any pan layer. The 
Cherenkov radiation angle in water is 35^, while it is greater in plexiglas (48^) and in 
glycerin (47^). 

The photon number in the wavelength interval (Ai = 3000A, A2 = 8OOOA) emitted 
in the radiator (n = 1.49, / = 5.5 cm) is given by the Frank and Tamm formula ^lOj : 

where a = 1/137; n is refraction coefficient; m^^E^ are muon mass and energy. 
Integrating it above the threshold energy {E^ > 0.25 GeV) with the muon spectrum 
measured at sea level we have = 2850. Then the ratio of the air Cherenkov light 
intensity Q{R) at distance R from the shower core to A^^ is 

Q{R)S ^ AsK; 

~ Ar{l-Rsy 

where Ag^Ar are the detector signal amplitudes from the shower and radiator; S is 
photocathode area; Rs < 0.04 is the reflectivity factor of air-glass junction. 

— _ !^:N,{X)S{\)KA\)dX _^^^^ 

j;:^ N,{\)s{\)d\ ■ ' 

where Kr{X) is a radiator transparency of effective thickness 4 cm (taking into account 
a radiation cone). The spectral characteristics of PMT, plexiglas and glycerin are 
measured with the spectrophotometer. Similar values were derived in the case of 
water radiator [9 . 

Optical radiators used in the laboratory for absolute calibration of the detector 
produce photon number insufficient for the routine ratio calibration of a multitude of 
array detectors in the field. Instead, the scintillator disk (D=15 cm, h=5 cm) is in 
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Figure 7. A setup to calibrate the air Cherenkov hght detector. 



use as the light source in this case. Intercahbration of different radiators has been 
performed using the experimental setup (figure [t]) . 

The total calibration error of 21% consists of uncertainties in: 

• number of photons at photocathode (10%); 

• measurement of signal amplitude from radiator (9%); 

• conversion factor from radiator to scintillator (7%); 

• signal variance throughout the photocathode area and acceptance angle (15%). 
5. Atmospheric extinction of light 

Cherenkov light undergo extinction in the atmosphere because of absorption on 
molecules by Rayleigh scattering and Mie scattering by aerosols and absorption. 
During the observation periods in winter (average temperature —40^ C) Cherenkov 
light absorption in the atmosphere is negligible and only molecular and aerosol 
scattering of photons are taken into account. 

Molecular scattering is almost constant, while aerosol concentration in the 
boundary layer above the surrounding terrain is of diurnal and seasonal variability. 
We are using the event rate of the showers in the energy range 10^^ — 10^^ eV detected 
with a PMT subset in hour and 15 min intervals in order to monitor the atmospheric 
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transparency for the light generated in EAS [3] . The method is based on the Cherenkov 
hght flux proportional to the primary energy: 

^ Ni{> Qthr) ^i/^ 

where 7V(> Qthr) is integral number of events with the light intensity above the 
threshold detected in z-th and basic periods; r is atmospheric extinction coefficient; n 
is the energy spectrum index. 

Resultant extinction coefficient is a product of variable due to current aerosol 
concentration and basic coefficient tq caused by the molecular scattering and minimal 
aerosol extinction. Rayleigh scattering parameters evaluated for the Yakutsk array 
conditions give = 0.9 [6 and the corresponding attenuation length A^^ = 300 km, 
averaged over wavelength range under consideration. 

The integral spectrum index of the Cherenkov light flux around the knee is [6]: 

K = 1.5 ±0.03, £;o < 3 X 10^^ 
n = 2.1 ±0.04, £;o > 3 X 10^^ 

Evaluating the coefficient tq as a function of Xmax in the shining point 
approximation [11 by assuming the attenuation length due to aerosol scattering 
constant above 1.5 km height and decreasing below [12 , one can derive the average 
extinction coefficient shown in figure [8] The approximation is 

51^30+ (1.11 ±^i2)x 10. - 

where Xmax is in g/cm^. 
Finally, we have found 

r = 0.765 ±0.015 

at ^ = 20^, £^ = 7 X 10^^ eV averaged over all observation periods, and inserting 
spectrum index in equation ([2|. 

6. Lateral distribution of air Cherenkov light 

During ~ 15000 hours of observation ~ 60000 showers of energy above 6 x 10^^ eV 
were detected by the medium Ci subset. The autonomous C2 array data consist of 
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Figure 9. Comparison of lateral distribution of air Cherenkov light. Triangles - 
this measurement; circles from ref 14 . 

~ 200000 showers with E > 1.2x 10^^ eV detected during ~ 3200 hours of observation. 
The Cherenkov trigger condition only is used to select these showers. 

Relative dispersion of the light intensity observed with different detectors at 
R 500 m from the shower core has been evaluated selecting showers with 3 detectors 
hit in the interval 400 < R < 500 m not far from each other (< 75) m. Intensity 
dispersion contains, in addition to inherent fluctuations instrumental errors, core 
location error, intercalibration error, etc. Supposing chance variation of the signal 
around an average Q{R) ex R~^-^^ [13], we have found SQ/Q = 0.25 ± 0.13. In order 
to minimize the Q{R) uncertainty due to EAS axis location error, one has to select 
showers with axes within R < Ropt m, (80 < Ropt < 250, if > 2 x 10^'^ eV; 
80 < Ropt < 400, if £^0 > 2 X 10^^ eV) where stations are spaced closely. 

In the highest energy domain measurements of Cherenkov light exist, carried 
out by the Haverah Park group J4] . In figure [9] our results are given in comparison 
with [14]. The data are normalized at 200 m core distance and exhibit consistency in 
the lateral distribution function (LDF) shape. The solid lines are an approximation 
to our data given in [9 . 

Model simulations give a variety of lateral distribution functions, some of them 
are close to our data, as illustrated in figure [lOj The lateral distribution depends 
on the attenuation lengths, multiplicity of secondaries in the interactions, primary 
mass composition etc., which can be parameterized with X^ax- Another considerable 
influence is the angular distribution of electrons in the shower. The combination of 
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Figure 10. Results of model calculation and the Yakutsk array data (circles). 
The lines illustrate results of 1: Hara et al. |7j, 2: Dyakonov [15| , 3: Ivanenko et 
al. [16], 4: Lagutin et al. |T7], and 5: Dedenko et al. [18]. The primary particle 
energy is indicated at the left of data/lines. 



these factors results in different LDFs. The results of Ivanenko et al. and Dyakonov 
are compatible with our measurements. 

Our cumulative results on LDF measurements (zenith angle < 30^) are given 



in figure [IT] The data of both subsets Ci , C2 are parameterized by the intensity at 
150 m from the shower core, Q150, the only core distance really present in the shifting 
range of measurements when the primary energy is rising from £^0 ^ 10^^ to 10^^ eV. 
The data are consistent with the previous results of the Yakutsk array concerning the 
Cherenkov light and can be described by the suitable EAS model simulation [151 [T6]. 

No abrupt change of LDF parameters is seen, so we choose a rather smooth 
approximation curve to fit the experimental data in the whole energy range [4^: 

where Ri = 60 m; R2 = 200 m; 6 = (1.14 ± 0.06) + (0.30 ± 0.02) x IgQiso- 
7. Energy estimation 

The total flux of Cherenkov light emitted, Qtot^ is our main estimator of the primary 
particle energy. In order to derive the relation between Qtot and ionization loss of 
the shower electrons in the atmosphere we have used formula ([T]) for the number of 
photons induced, A/"^, in ratio to that from optical radiator, N^: 

dN^ dNr 1 / 1 \2\ 

~dh ~ ~dhl-l/nl^ 



where v is electron velocity; n, are the refraction indexes in air and radiator, 
respectively; c = 1. 
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Figure 11. Air Cherenkov light radial distribution. Our experimental data are 
given by the black and white points alternately, not to confuse the adjacent data. 
The curves are approximations according to Equation with Q150 fitted to the 
particular data. An axis distance R = 150 m is indicated by the dotted line. 



Emission is possible when vn > 1 or E'^ > E^^^ = m^/(l — 1/n^). The number 
of photons induced in the depth interval dx is 

5^ = C(i-§^), (5) 

where ( = ^ 1-^1/ n"^ '> Pr is radiator density; po^^o are air density and refraction 
at the observation level. 

The total light flux is given using the electron differential energy spectrum in the 
shower at the depth x: 



{xo) = / dxT{xo - x) / dECi^ 

Jo JEthr ^ 



dN{x,E,Eo) 
' dE ' ^ ' 

where r{xo — x) is the extinction coefficient of light along the path xq — x. 

The number of electrons in a shower is approximately N{x^ E^ Eq) 
N{x,Eo)x{E), where N{x,Eo) is the number of electrons with ^ > 0; x(^) is the 
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Table 1. The fraction of air Cherenkov hght, Ag, charged particles, As, and 
muons, A^, actually measured in EAS events in the Yakutsk array experiment. 
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universal electron spectrum at the shower maximum. Then 

Qtot = 2CE^hr / dxr{xo - x)N{x, Eo) / x{E)-^ (7) 

Jo JEthr ^ 

leads to the relation with the ionization loss of electrons in the atmosphere Ei = 
£0 j^' N{x,Eo)dx/tQ: 

Ql=^ ^ 2^0(0.275 - 0.283 ^^^"^^"^""^ )^, (8) 

^0 ^0 

if we assume the extinction r = 1 [T9l. 

In the real case of the Yakutsk array conditions T = -30^ C, P = 754 Torr and 
the extinction of light described in Section 5, it was found that the relation is [19 : 

E^/Qtot ^ (3.01 ± 0.36) X 10^(1 - X^,,/(1700 ± 270)) (9) 

for Xmax e (500, 1000) g/cm^; the accuracy of the relation is ~ 5%. It demonstrates 
the advantage of the air Cherenkov light measurement technique: the relation is 
determined by Xmax and extinction of light only; interaction model dependence is 
parameterized by means of X^ax- 

In [Appendix A we analyze the shower parameters governing the energy fractions 
transferred to EAS components. 



7.1. Experimental evaluation of the energy transferred to EAS components 

The energy fractions of the main EAS components can be estimated using the Yakutsk 
data. The ionization loss of electrons is measured when detecting the total flux of the 
Cherenkov light at ground level. The detector disposition of the Yakutsk array is 
appropriate to measure Qtot in the range above about Qi50 = 10^ m~^ as is shown in 
table [1] In each Q150 bin the LDF extrapolation formula Q is used to calculate the 
total flux. 

Conversion of the measured Q^^^ to Ei is carried out along equation (|9| with 
parameters relevant to the particular observation period, taking into account detector 
calibration and atmospheric extinction of light. 

Other portions of the energy carried out by electromagnetic and muonic 
components beyond sea level is evaluated via the total number of electrons: 

Eg = CoNeXe/toj 

where eo, to are the critical energy and radiation length of electrons in air; attenuation 
length Ae is derived from zenith angle dependence of A^e [20^ ; and muons measured at 
the ground level: 



Measuring EAS with Cherenkov light detectors 



13 



Table 2. The primary energy fractions of EAS components: ionization loss of 
electrons in the atmosphere, Ei, and in the ground, Eg] energy of muons at the 
ground level, E^] energy of EAS components unobservable at the Yakutsk array, 

Eunobs- Eq = Ei -\- Eg + Ell + Eunobs- = 0^. 
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Figure 12. The electromagnetic component energy estimation from the Yakutsk 
array data: Ci subset (squares), C2 subset (rhombuses); solid curve is the result 
of CORSIKA/QGSJET based estimation. 



where the average energy of muons E^ is taken from the MSU array data [2T]. 

Residuary energy fractions transferred to neutrinos Ei^^ nucleons Eh etc., 
unmeasurable with this array, are estimated using computational modeling [19 . The 
resulting apportioning of the primary energy £^0 is given in table [2] 

The energy fraction carried by electromagnetic component (£^em = Ei ^ Eg) 
appears to be the basic contribution to the total energy of the shower, and its 
energy dependence (measured with Cherenkov light detectors + scintillators of the 



Yakutsk array [22]) is illustrated in figure [12] in comparison with CORSIKA/QGSJET 
estimation [23^. 

Due to the air Cherenkov total light fiux and the electron and muon number which 
are experimental values measured at ground level, only about 5% of the primary energy 
in the interval £^0 ^ (10^^ 7 10^^) eV is calculated by using model assumptions. So we 
consider the energy estimation used to be model-independent within these bounds. 

Moonless nights, when air Cherenkov light measurements are possible, constitute 
~ 10% of the observation period. In order to evaluate the primary energy of the bulk 
of showers, the correlation 

^600 = 1.56 X lO-^Qls^o' 



is used between the charged particle density at 600 m from the shower core. Si 



600, 



and the light intensity at 150 m from the core (figure 13) which, in turn, is related 
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Table 3. Parameters of the relation between Q150 and the total flux of air 
Cherenkov light and the energy of the primary particle initiating EAS. 





ai 
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< 10^ 
10^ - 10^ 
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>109 


2.39 X 10^ 
1.33 X 10^ 
7.91 X 10^ 
5.75 X 10^ 


0.90 
0.94 
0.97 
0.98 


(6.87 ± 1.44) X 10^^ 
(3.78 ±0.72) X 10^0 
(2.90 ±0.29) X 10^0 
(5.59 ±0.56) X 10^ 


0.87 ±0.02 
0.91 ±0.02 
0.93 ±0.02 
1.01 ±0.02 
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Figure 13. The correlation between charged particle density at 600 m from the 
core, Sqoo, and air Cherenkov light intensity at 150 m, Q150, measured in the 
same showers with < IS*-*. 



to the total flux of the Cherenkov light in the atmosphere Qtot — <^iQi50 5 parameters 
are given in first columns of the table [3j Using derived relation between the primary 
energy and Qtot (equation [9| the final formula was found to link Qi50 with £^0 in the 
interval 6> < 15° [l^: 

^0 = a2Q'ilo, (10) 

where the numerical values of a2, 62 are summarized in the last two columns of table [3] 
for different Qiso intervals. 

The observed densities S'soo/^'eoo at various zenith angles are connected to the 
'vertical' one {Q = 0°) along attenuation curve [24]. In order to measure the 
attenuation length of these densities for fixed energy, we have used two different 
methods - well-known equi-intensity cut method, and fixing the Cherenkov light 
intensity at 400 m from the core as the equivalent of the primary energy, taking 
into account the light absorption in the atmosphere. In figure 14 the results are given. 
Experimental points are consistent with each other for the two methods used and can 
be described by the sum of two components - a soft component (electrons, attenuation 
length Ae = 200 g/cm?) and a hard component (muons, = 1000 g/cm?) [20 : 

S^m{0) = 53oo(0°)((l - /330o)exp(^^) + p^^^exp^^^)), (11) 

where /^aoo is the hard-component fraction. Attenuation curve for ^eoo is the same 
but (5 is different: 

/?3oo = (0.563 ±0.032)53oo(00)-°-i85±o.02^ 
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Table 4. The highest energy EAS events detected with the Yakutsk array 



Date 


0' 




6Eo,% 






18.02.04 
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140.2 


07.05.89 


58.7 


20.14 
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161.6 


21.12.77 


46.0 


20.01 


40 


50.0 


220.6 


15.02.78 


9.6 


19.99 


32 


15.5 


102.0 



f3eoo = (0.62 ± 0.006)56oo(0°)-° ™ °^ (12) 

Experimental uncertainties when estimating the EAS component energies are 
summarized in the last column of the table |2j The main contribution arise from SEi 
which is governed by uncertainties in the atmospheric transparency (15%), detector 
calibration (21%) and the total light flux measurement (15%). Errors in estimation of 
A/'e, Ae, A/"^ determine the next two items (for ionization loss in the ground and SE^). 
Resultant uncertainty in energy estimation is the sum of all errors weighed with the 
second/third columns of the table [2j SEq ~ 32% [19 . Extra 20% are added due to a 
Qi50-5'6oo conversion uncertainty. 

To illustrate the energy estimation method used, four showers detected at the 
Yakutsk array in the range Eq > 10^^ eV, < 60^, and axes within the array area, are 
given in table |4] (one event is added slightly below the threshold because the energy 
estimation error is larger than the tiny difference). 



8. The energy spectrum of cosmic rays derived from air Cherenkov hght 
measurements 

The Cherenkov light detector subsets of the Yakutsk array give us the opportunity to 
reconstruct the energy spectrum of cosmic rays in the energy range from Eq ~ 10^^ to 
6 X 10^^ eV. The total number of electrons and muons measured at the ground level 
are used to estimate the additional energy fractions carried by EAS components but 



the final relation (10) comprises the light intensity alone. 
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Figure 15. Acceptance area of air Cherenkov light detector subsets of the 
Yakutsk array. The sohd curve is for autonomous, C2, and the dotted one for 
the medium, Ci, sub-array area. 



The intensity of CRs has been evaluated using the number of EAS events derived 
from the showers detected with PMTs, and aperture SeffTVt of the sub-arrays in 
a particular energy interval, where the array area bounded by the perimeter, Seff^ 
depends on the primary energy and zenith angle; Vt is the acceptance solid angle; and 
T is the sum of observation periods. Under the condition of the array configuration 
changing due to rearrangements and non-active detectors, Seff has to be calculated 
for a given period using the Monte Carlo technique. 

The acceptance area has been simulated as a function of Qiso averaged in the 
zenith angle interval (0^,30^) for the two Cherenkov light detector subsets shown 
in figure [l] A lateral distribution fit accord to Eq. Q is used together with 
instrumental and statistical errors (Gaussian with 25% relative deviation) to model 
the trigger of each detector subset with 100000 fake showers. In the case of the 
autonomous sub-array, the Cherenkov trigger is simulated, while in the Ci case the 
double trigger for scintillator and PMT signals has been modelled. The average 
{Seff = Sontriggered/^^'^i whcrc 5*0 is the array area inside the perimeter; ritriggered 
is triggered number of events) is shown in Figure 15 versus the parameter Qiso- 
Corrections for inoperative detectors are not shown here. 

The shower data gathered after the latest array re-configuration were used to 
work out the spectrum. Namely, 1993-2007 for the medium Ci subset, and 1995-2007 
for autonomous subset. In order to evaluate the intensity of the primary flux we have 
collected data during observation periods with a light extinction better than 0.65 and 
shower axes within area of the corresponding sub-array. 

The resulting differential all-particle spectrum of cosmic rays is shown in 
figure 16 in comparison with the data from other Cherenkov detector arrays, namely. 



BLANK A 1 25 and Tunka [26]. The present data (given in tabular form in Appendix 
[c{ tables CI, C2) exhibit the spectrum irregularity near E ~ 10^^ eV, the 'ankle'. 



seen by all arrays in the area [27J; at lower energies the results of all three arrays 
are compatible with a 'knee' at £^ ~ 3 x 10^^ eV revealed in the pioneering works 
of the MSU array 1281 . Below E = 10^^ eV there is a transition region between the 
two subsets, Ci,C2, so the 'second knee' visible here may be due to the data sewed 
together. 

Hereinafter, two energy regions below and above E = 10^^ eV are analyzed 
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Figure 16. Cosmic ray flux measured by arrays equipped with air Cherenkov 
hght detectors. Vertical bars indicate statistical errors. 



separately in order to consider the two irregularities in the spectrum. 

The region below 10^^ eV is shown in detail in figure 17 where CR intensity is 



multiplied by E'^-^^ in order to emphasize the spectrum irregularity. In addition to 
spectra around the knee measured by Cherenkov detector arrays (BLANKA, Tunka 
and Yakutsk taken from the previous Figure), the data of scintillation detector arrays 
(Akeno [29 , KASCADE [30^ and Tibet [ST) are given for comparison. Due to air 
Cherenkov light signal proportional to the number of electrons in the shower, no 
difference is expected in the spectrum shape measured with scintillator or PMT arrays. 
Measured CR intensities are corrected in the case of Cherenkov detector arrays along 
the algorithm in [Appendix B| Integral energy spectrum indices below and above 
the knee are assumed to be 1.67 and 2.1, respectively, as measured by the Tibet- 
Ill array, while energy evaluation errors are estimated as 0.12 (Blanka), 0.2 (Tunka) 
and 0.25 (Yakutsk). Original intensities from Akeno, KASCADE and Tibet arrays 
are not changed because the energy /intensity reconstruction procedures from N^^N^ 
measurements include the conversion factor ( |B.1[ ) needed. 

For the reconstructed energy and intensity there is an interaction model/primary 
composition dependence in all the data from arrays. For instance, the Tibet-Ill results 
indicate a ~ 20% systematic error due to chemical composition and ~ 10% discrepancy 
between QGSJETOlc and SIBYLL2.1 interaction models below E = 10^^ eV. The 
uncertainty may by much worse if to use interaction models not carefully tuned to the 
measured EAS observables. 

We have used energy correction factors to compare the measured differential 
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Figure 17. The energy spectrum of cosmic rays around the knee region. 
Observational data are from arrays equipped with Cherenkov hght and 
scintihation detectors. Two panels show the spectra before (left) and after (right) 
energy corrections applied to the data, as described in the text. 



spectra around the knee. The intensities (J x E'^-^^) shown in the left panel of the 
figure [17] were shifted by factors 1.2, 1.05 and 1.05 in Blanka, Tunka and Yakutsk 
cases, respectively; other three results were shifted along energy correction factors 
0.95 (Akeno), 0.97 (KASCADE) and 1.02 (Tibet). The resultant spectra are plotted 
in the right panel of the figure [T7[ All data exhibit the knee approximately in the 
same energy interval while the intensity /energy estimates are different especially in 
the Akeno and BLANCA fiux data. 

As examples of model calculations, four predicted spectra are shown together 
with the Cherenkov light data of the Yakutsk array in figure [18] Only the shape of 
the all-particle spectra can be compared with the data because of the free parameters 
in models - the intensity of CRs and, to a lesser degree, a knee position due to the 
magnetic field uncertainty in the sources and in the interstellar medium. 

All the model spectra are compatible with experiment, especially, in view of 
the dispersion in the results of arrays (e.g. figure 17) which is greater then the 
shape variation due to models. At energy £^ > 0.1 EeV the anomalous diffusion 
model spectrum is harder then that observed in Yakutsk, but the points are within 
experimental errors. The lack of CRs in SNR acceleration and diffusion models above 
energy 0.1 EeV can be filled up by the extragalactic component. 

There is seen to be a hint of the fine structure in the energy spectrum measured 
with the Yakutsk array. Although it is possible that the undulations around the knee 
are caused by the instrumental error^ but on the other hand, there are 40 size spectra 
and 5 Cherenkov light spectra measured before 2005, which demonstrate the second 
excess ('peak') at lg£^ = IgE^^nee + 0-6 besides the knee itself [36 . The second peak 
in our data is approximately at this energy. However, more Cherenkov light data are 
needed to scrutinize the subject. Our autonomous sub- array is able to supply with a 
sufficient sample of EAS events within the next few years. 

Recently, Lagutin et al. [37J examined the contribution of a nearby SNR-type 
source to the energy spectrum of CRs produced in anomalous diffusion model. They 
found several sequential peaks caused by H, He and CNO nuclei around the knee in 
all-particle spectrum, confirming results of the single source model [34] in the case of 
the 'background' anomalous diffusion model. 



X vertical bars show statistical errors only 
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Figure 18. Our energy spectrum and the results of galactic cosmic ray 
simulations. Curves: SNR acceleration model [32] (BV); diffusion model [33] 
(KP); single source model |34] (EW) and anomalous diffusion model 35 (LNU). 



Table 5. Energy estimation errors, cr, and intensity conversion factors, Rj, for 
EAS arrays. 



Array AGASA HiRes PAO Yakutsk 

25 17 22 32 

Rj{lgE < 18.5) 1.18 1.08 1.14 1.31 

i?j(18.5 < IgE < 19.8) 1.12 1.05 1.09 1.20 

RjjlgE > 19.8) 1.65 1.26 1.47 2.27 



Our conclusion concerning the part of the energy spectrum below 10^^ eV is that 
while all four models considered are compatible with our measurements in the knee 
region, e.g. the intensities below and above the knee, the only model able to describe 
the fine structure in the spectrum is the single source model. So this model combining 
a recent nearby SNR with the background (anomalous) diffusion of CRs in Galaxy is 
the best fit for the Yakutsk array data. 

When comparing the upper half of the spectrum measured for E > 10^^ eV, with 
the energy spectra observed by giant EAS arrays, including that of the Yakutsk array, 
the somewhat different energy estimation has to be taken into account [3l [20l [38] . 
Due to the comparatively small acceptance area of the Cherenkov detector subsets 
(figure 15), our data are reliable up to ~ 10^^ eV; this includes the ankle region. 
Another purpose is to compare UHECR intensities measured with different techniques. 

In the table [5] are given the energy estimation errors of the AGASA [39] . 
HiRes [40] , PAO [41] and Yakutsk [19] experiments and the intensity conversion factors 
calculated as described in [Appendix B| for energy bins where the index is constant and 
the variation of instrumental errors is negligible. The integral energy spectrum index 
of CRs is assumed to be 2.3 and 1.9 below and above the ankle in the spectrum, and 
Ai: = 4 at Ig^ > 19.8, along the results of HiRes fitted by the broken power law [40] . 

Observed UHECR spectra are given in the left panel of the figure [19] with the 
intensity correction factors applied: observed intensities are decreased by Rj and 
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Figure 19. UHECR energy spectra measured by AGASA (open rhombuses) |39| , 
HiRes I, II (open and filled triangles) PAO (filled rhombuses) gl], the 

Yakutsk array scintillation detectors (open, crossed and filled circles) [20], and 
air Cherenkov light detectors (squares) 43 . Statistical errors are shown by the 
vertical bars, while horizontal bars indicate energy bins. The spectra are given 
before (left panel) and after (right panel) energy corrections applied to the data. 



Table 6. Correction factors to energy sea 
Re, averaged in the region E > 10-*^^ eV 



factors to energy scales for the pairs of EAS arrays/detectors, 
region E > 10^^ eV. 





AGASA 


HiRes 


PAO 


Ysc 


Ych 


AGASA 


1 


0.75 


0.63 


1.05 


0.82 


HiRes 


1.33 


1 


0.85 


1.40 


1.08 


PAO 


1.6 


1.2 


1 


1.70 


1.30 


Ysc 


0.91 


0.71 


0.6 


1 


0.75 


Ych 


1.22 


0.93 


0.80 


1.33 


1 



spectra are displaced over 0.434(7 along lg£^j§] This is preliminary, crude procedure to 
reconstruct the spectrum breaks, but for the purpose of the intensity comparison it 
may be sufficient. 

Two spectra of the Yakutsk array are compatible within errors above 10^^ eV 
but diverge at lower energies. The possible reason of a discrepancy can be systematic 
errors in primary energy and CR intensity estimation near the scintillator threshold 
of the trigger-500. The work is in progress to surmount the divergence. 

Comparison of UHECR spectra measured using different detectors and energy 
reconstruction methods infer the existence of systematic differences between resultant 

energies obtained at the arrays. The average values, Ei^ should be corrected, too||i] 
While the true primary energy is unknown, cross calibration of energy estimation 
methods can be carried out adjusting correction factors, Re^ for the pairs of observed 
spectra, converging it together. Then the resulting spread of factors elucidates the 
confidence interval for the CR energy estimated. 

Table [6] demonstrates the variety of correction factors for the energy estimation 
methods of the different groups. 

To illustrate the result of corrections Rj x Re applied to estimated energies and 

§ This may be unwarranted in the PAO case due to (42]. 
II for the given primary energy E 
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Figure 20. The Yakutsk array data from the previous figure (with energy 
correction) in comparison with the results of UHECR propagation models given 
by curves: Berezinsky, Gazizov & Grigorieva (BGG); Bahcall & Waxman (BW); 
Wibig &; Wolfendale (WgW); De Marco Sz Stanev (MS). 

intensities, the measured spectra are re-plotted in the right panel of figure [19] Energy 
scale factors are used here (arbitrarily, but preferably close to the median) from the 
last column of table [6j although any other column may be used as well. 

The shape of the UHECR spectrum measured by all the arrays is compatible 
within errors, if the energy estimations are calibrated. Namely, the observed position 
of the ankle and energy threshold of GZK suppressiorjj] are in satisfactory agreement. 

A difference in the energy of EAS primary particle estimated basing on the data of 
the two subsets of the Yakutsk array detectors - scintillators, Ysc^ and air Cherenkov 
light detectors, Ych^ originates, presumably, in the systematic error of the relation 
used between the charged particle density at 600 m from the shower core and the light 
intensity at 150 m from the core. 

The UHECR propagation modeling results are illustrated by the four examples 
in figure [2Q| in comparison with the Yakutsk array data. 

Bahcall and Waxman gave two-component (Galactic+extra-galactic) model with 
the shape insensitive to the choice of absolute energy scale [44]. 

Another approach was used by Berezinsky et al. assuming UHECRs as 
extragalactic protons from uniformly distributed sources [45 . The electron-positron 
pair production in collisions of protons with relic photons results in the energy 
spectrum of extragalactic CRs with the 'dip' feature in this model. 

De Marko and Stanev' model fits the UHECR spectra measured by AGASA 
and HiRes with different injection spectra at CR sources that are uniformly and 
homogeneously distributed in the Universe [46 . The best fit they found assuming that 
cosmic rays {E > 10^^ eV) are protons, varying the index, emissivity and cosmological 
evolution parameters of the injection spectrum, is given in the figure. 

Wibig and Wolfendale focus on the ankle in the primary energy spectrum 
attributing it to the rapid transition from Galactic to Extragalactic component of 
CRs [47]. The sum of a smoothly falling Galactic spectrum and a power- low EG 
spectrunj^ fitted to the Yakutsk array scintillator data is shown. 

^ except AGASA data 
+ index -2.37 



BGG 

BW 



WgW 

MS 
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All the models used demonstrate GZK suppressiorj^ and the ankle features in 
agreement with the data from arrays. It is not surprising in view of the fact that the 
source emissivity and the ankle position in the energy scale are free parameters, and 
GZK effect is embedded in models. 

The lack of EAS events above 10^^ eV observed with Cherenkov light detectors 
hinders in choosing a model of the better fit. The index of the energy spectrum 
observed in the interval (10^^,3 x 10^^) eV is d\gJ/d\gE = —2.1 for Cherenkov light 
detectors data and is —2.7 for charged particle detectors data. If we use these values 
as the confidence bounds then model indices are within the interval. The model of 
Bahcall and Waxman is closest to our Cherenkov light detectors data and can be 
considered as preferable among equally matched. 

9. Conclusions 

The total flux of air Cherenkov light with subsidiary data on the electron and muon 
sizes at the ground level is used to estimate the energy of the primary CR particle 
initiating EAS. The relation between the total flux and ionization loss of electrons 
in the atmosphere is derived which depends on the extinction of light and X^ax 
parameters; the latter is the only parameter to accumulate the interaction model 
dependence of EAS development in this case. 

The independent measurement technique based on the Cherenkov light detectors 
of the Yakutsk array enabled us to observe the cosmic ray energy spectrum in the range 
from E ~ 10^^ eV to 6 x 10^^ eV. Two spectra measured with different detectors of the 
Yakutsk array - scintillators and Cherenkov light detectors, exhibit an ankle feature 
below E = 10^^ eV. The autonomous sub-array data confirm the previous observations 
of the knee at £^ ~ 3 x 10^^ eV. A comparison of our results with the data of other 
EAS arrays shows the compatibility of spectra if the energy estimations are corrected. 

The energy spectra predicted for several models of galactic CRs demonstrate 
agreement within experimental errors with our data in the knee region. However, 
only the single source model describes the fine structure of the spectrum observed 
attributing it to the contribution of nuclei from a recent nearby supernova. We have 
chosen this model as the best fit for our data below 10^^ eV. 

At the highest energies extragalactic CRs forming a dip, or transition between 
galactic and extragalactic components above 10^^ eV are thought to be responsible for 
an ankle detected in the energy spectrum. A comparison of our data with models of 
these types shows a satisfactory agreement with all of them (the best fit with Bahcall & 
Waxman's model), so we cannot distinguish between different scenarios of the ankle 
formation basing on the energy spectrum measurement alone. There is a need for 
additional data, presumably, concerning the mass composition of UHECRs, in order 
to elucidate the origin of an ankle in the energy spectrum. 
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Figure Al. The energy carried in a cascade by nucleons (A^), charged pions (tt), 
muons and neutrinos (fi -\- u), electrons and photons (e + 7). Dashed curves are 
analytic expressions with Xi,K]y = const, B-jr = 0: Ejy = Eoexp{—K]\jx/XN); 
E^ = 2/3£;o(l - exp(-KNx/XN))exp(-x/X^/3y, Ee^ = Eq - En - E^. 
Solid curves are (5-model results: cross sections are supposed rising oc 0.08^n£'; 

= 120 GeV; = 0.5; Us oc E^/^; Eq = 10^^ eV. The open circle represents 
the asymptotic estimation for Ee^ of the CORSIKA/QGSJET code. 
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Appendix A. Energy balance of EAS components 

The energy fractions of the EAS primary particle transferred to the shower components 
can be described on the basis of hadron transport equations. If Ej.^ {k = 
TT, iiv^ 67) is the energy transferred to nucleons, charged pions, muons+neutrinos, 
electrons+photons, a few cascade parameters determine the ratios between E^ - i.e. 
the energy balance in the shower. For instance, the transport equation for the charged 
pions density 7r{x,E) at depth x is: 

d7r{x,E) ,1 ^TTw Z7^ ^ 



dx 

+ ^ / N{x,U)w^N{E,U)dU, (A.l) 

where the interaction mean free paths A7^,AAr are assumed to be constant; 
Wt^t^{E^U)^Wt^n{E^U) are the spectra of charged pions produced in pion-air and 
nucleon-air interactions, can be transformed (integrating over E'^) to: 

dE^ ^ _E^ _ BMx,E>0) ^ 2^ ^ 2KnEn 
dx X SAtt SAat 

where E^{x) = J^"" 7r{x, E)EdE; 7r{x,E > 0) = J^"" 7r{x, E)dE; Kn is the nucleon 
inelasticity assumed to be constant; En = Eoexp{—KNx/XN)' 

In the energy range E ^ B^^ the only parameters to define the solution are 
Kn/Xn and A^^. It means that the energy transferred to charged pions is independent 
of the spectra of pions produced in nuclear interactions. Hence, in the general case, 
we can use simple (5- model with the production spectrum Wik{E^ U) — ns5{E — U /ris)^ 
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where is the multiphcity of secondaries, to balance the components energy in a 
shower. Because of the net value of E^j^y/EQ < 0.1, the uncertainty due to the 
simplified model should be of second order of magnitude. 

To summarize, the main model parameters governing the energy balance in 
the shower are the average inelasticity coefficients, mean free paths, multiplicity of 
secondaries and the fragmentation rate of the primary nucleus. The influence of other 
model characteristics such as 'the form of the rapidity distribution of constituent 
quarks' in collisions is weak. 

The exact analytic solution of equations in the case of constant X^^, Xn, Kn 



and = is shown in figure Al together with (5-model results. The letter is 
demonstrating the influence of rising cross-sections and the real decay rate of charged 
pions. An asymptotic {x = oo) estimation of E^^ with the CORSIKA/QGSJET code 
at Eq = 10^^ eV [23 is shown in the figure as well. 

Appendix B. Conversion of the measured cosmic ray intensity to the 
spectrum of the primary beam 

There is a correction to be applied to the measured intensity of CRs before any 
comparison of the energy spectra observed by different EAS arrays. 

The quantity E = 'primary particle energy' that has been estimated after a shower 
detection, and the actual energy of the particle, E, which has initiated the EAS, are 
different values, connected with each other by a relation to be found. Estimated energy 
has distribution around the given mean value g{E^E)^ formed by the instrumental 
errors and fluctuations of the shower parameters with a RMS deviation, a. Energy 
fluctuation is small in comparison with instrumental errors. Our aim here is to 
calculate the exact difference between the observed intensity of cosmic rays, J{E)dE, 
and the original one J{E)dE in the case of a rapidly falling power law spectrum. 

The problem was solved, in general, by Zatsepin and Kalmykov in the previous 
century. The measured number of EAS particles, the so-called shower size, A/'g, is 
connected with the primary energy. The function g{Ne,E) depending on E/N^ has 
been found by Zatsepin [48 . Then Kalmykov has calculated the measured intensity 
in the case of a lognormal distribution of [49 : 



J{Ne) = MNe)exp{ \ , (B.l) 



where ctat is RMS deviation of In A/'g; is spectrum index; ajy = -^^j^- 

Here we also assume the lognormal distribution ofy = \nE with the average value 
equal to \nE. The observed intensity of cosmic rays is then given by the convolution 
of the primary spectrum, J{z) = Jo exp(— and the distribution of instrumental 
errors and fluctuations 

poo ^rj^rg(_y^\ 

J{z-y)g{y)dy = Jo / exp{-f^z^f^y)^^)^^dy.{B.2) 
-oo J-oo y Zira 

The resultant initial-to-observed intensity conversion factor is 

2 2 

Rj = J{z)/J{z) = exp{^). 



The necessary conditions are a constant index and RMS error, or at least both changing 
only slowly with energy. 
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Figure Bl. Simulation results of the spectrum measurement with ankle and 
knee features. 



A distinct feature of the break in the spectrum is its shift along the energy scale. 
Due to the smeared transition between the conversion factors below and above the 
break, its position moves upward in energy. 

To illustrate this we have simulated the primary energy spectrum with ankle and 
knee: random primary energies have been generated according to a broken power law 
with the index 7 = 3.3 and 2.9 below and above 10^^ eV in the case of ankle, and 
7 = 2.9 and 5 for the knee. The energy estimation is modelled by adding randomly a 
Gaussian error to \iiE with a = 0.32. Figure [Bl] shows the 'measured' and primary 
spectra. The different shifts in the logarithm of intensity on both sides of the break 
are clearly seen, as well as that the position of the observed break moves to the right. 
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Table CI. Cosmic ray flux measured with the C2 sub-array detectors. Statistical 
errors are given where appreciable. 

Energy, eV J(£^), sr~-^ s~-^ eV~-^ 
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Table C2. Cosmic ray flux measured with the Ci sub-array detectors. 
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EeV 


J{E), m-^sr-^s 
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